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ABSTRACT 

, Perturbation theory makes it possible to calculate the probability distribution function 

• (PDF) of the large scale density field in the small variance limit, a <C 1. For top hat 

smoothing and scale-free Gaussian initial fluctuations, the result depends only on the 
linear variance, ou near , and its logarithmic derivative with respect to the filtering scale 
-("linear + 3) = d log of inear /d log £ (Bernardeau 1994a). 
C*~) , In this paper, we measure the PDF and its low-order moments in scale-free simu- 

■ lations evolved well into the nonlinear regime and compare the results with the above 

predictions, assuming that the spectral index and the variance are adjustable parame- 
ters, n er r and a er r = u, where a is the true, nonlinear variance. With these additional 
degrees of freedom, results from perturbation theory provide a good fit of the PDFs, 
even in the highly nonlinear regime. The value of n g is of course equal to ni; near when 
tr«l, and it decreases with increasing a. A nearly flat plateau is reached when a 3> 1. 
i-^ ' In this regime, the difference between n e e and rqi near increases when nii noar decreases. 

For initial power-spectra with n.vt nea r = —2, —1, 0, +1, we find n c g ~ —9, —3, —1, —0.5 
when er 2 ~ 100. 

5h ' It is worth noting that — (3+n er r) is different from the logarithmic derivative of the 

r/3 , nonlinear variance with respect to the filtering scale. Consequently, it is not straight- 

forward to determine the nonlinearly evolved PDF from arbitrary (scale-dependent) 
initial conditions, such as Cold Dark Matter, although we propose a simple method 
that makes this feasible. 

' Thus, estimates of the variance (using, for example, the prescription proposed by 

J-j | Hamilton et al. 1991) and of n e s as functions of scale for a given power spectrum makes 

it possible to calculate the local density PDF at any time from the initial conditions. 

Key words: cosmology: theory - galaxies: clustering - methods: numerical - meth- 
ods: statistical 
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1 INTRODUCTION 

It is generally believed that the large-scale structures of the 
Universe, such as those seen in the distribution of galax- 
ies, arose through gravitational instability from primordial 
seed fluctuations. On scales larger than roug hly 8 /i _1 Mpc£| 
the dynamics of the Universe is dominated by gravitational 
forces, implying that the matter behaves as if it were purely 
nonbaryonic and collision-less. Thus, when curvature can be 
neglected (for scales below 500/i _1 Mpc), the evolution of 
the Universe can be described by the Vlasov-Poisson sys- 

* We adopt h = H /(100 km/s/Mpc), where H is the Hubble 
constant. 



tern of equations, in coordinates comoving with the cosmic 
expansion. In general, nonlinear coupling makes it virtually 
impossible to solve these equations analytically. Thus, on 
small scales, where density contrasts are large, iV-body sim- 
ulations are typically used to obtain numerical solutions. 
However, on scales larger than about 10ft _1 Mpc, the ampli- 
tude of the density fluctuations is modest and Perturbation 
Theory (PT) can be used to compute statistical properties 
of the cosmic fields. 

Here, we focus on the Probability Distribution Function 
(PDF) of the local density smoothed with a spherical filter 
of radius I at a time t, P(p,£,t). The PDF, a simple to 
measure statistic, has been widely used to characterize the 
statistical properties of the large scale galaxy distribution 
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(see, e.g., Alimi, Blanchard & Schaeffer 1990; Maurogordato, 
Schaeffer & da Costa 1992; Szapudi, Szalay & Boschan 1992; 
Bouchet et al. 1993; Gaztanaga 1992, 1994; Szapudi, Meiksin 
& Nichol 1996). 

Our aim is to determine whether there is a simple rela- 
tionship between the PDF obtained for the matter distribu- 
tion in the nonlinear regime and the one predicted by PT. 
The latter should, in principle, be valid only in the weakly 
nonlinear regime. In our calculations, we assume i) that the 
Universe is flat with Q = 1, and ii) that the initial fluctu- 
ations of the density field were Gaussian, pressure-less (i.e. 
cold), and scale-free. 

In real galaxy catalogs, the interpretation of the results 
is complicated by systematic effects, such as the so-called 
bias between the galaxy distribution and the matter distri- 
bution (e.g., Fry & Gaztanaga 1993, Juszkiewicz et al. 1995), 
redshift distortions in three-dimensional catalogs (e.g., Mat- 
subara & Suto 1994; Hivon et al. 1995) , and effects of projec- 
tion along the line of sight in two-dimensional catalogs (e.g., 
Bernardeau 1995). In what follows, we focus on the statis- 
tical properties of the matter distribution and thus neglect 
these effects. 

Let a 2 (£) be the variance of the PDF. In the weakly 
nonlinear regime (a 2 <C 1) the scale-free nature of the initial 
conditions implies 
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In general, the PDF of the local density can be characterized 
by the behavior of its cumulants (i.e. the connected part of 
the moments), {5 q ) c (see e.g. Bernardeau 1994a, hereafter 
B94 for a more precise definition). Note that {S 2 ) c — {5 2 } — 
a 2 (£), where 3 is the density contrast S = p/{p) — l. When the 
variance is small, the dominant contribution to all the cumu- 
lants can be calculated analytically from PT for a top-hat 
smoothing window (B94). The result is that the quantities 
S q defined by 



T 2(9-l) • 



(2) 



are all finite constants in the limit of small variance^]. The 
next-to-leading order term is expected to be proportional to 
a . so that 



S q {o) 



+ 0(a 2 ). 



(3) 



Comparisons with numerical simulations have shown that 
the PT predictions are accurate as long as a < 1 
(Juszkiewicz et al. 1993, 1995; B94; Baugh, Gaztanaga & 
Efstathiou 1995; Gaztanaga & Baugh 1995; Lokas et al. 
1995a,b). In this regime, the whole PDF can then be cal- 
culated analytically (B94), as a function of £, <rii ncar and 

^linear ■ 



Bouchet et al. 1993, Gaztanaga 1994, Szapudi, Meiksin & 
Nichol 1996), given possible statistical uncertainties (e.g., 
Szapudi & Colombi 1996) and other systematic effects dis- 
cussed earlier in this section. 

For power-law spectra, the set of Sq T can be deter- 
mined as a function of nii ncar (Juszkiewicz & Bouchet 1992; 
Juszkiewicz, Bouchet & Colombi 1993; B94; Bernardeau 
1994b): 
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and so on. 

In what follows, we use these definitions as fitting func- 
tions of the S q ratios, leaving nii ne ar as a free parameter. One 
can then invert the above equations to deduce, for each q, 
the value of niincar which yields the measured value, 5'™ cas , of 
the parameters S q . We use this approach to define functions 
n q (a), (q > 3), by 



Sq T [n q (a)] 



(a). 



(8) 



Of course, as discussed above, the measured n ? 's are ex- 
pected to have identical values in the weakly nonlinear 
regime (a <C 1), and be equal to the value of the initial 
index ni in car. 

In the nonlinear regime (a J; 1), however, results 
from perturbation theory are no longer valid, and we do 
not expect the functions n q to take a universal (constant) 
value. Actually, there is no a priori reason why the relation 
n q (a) = n q i{a) should hold for q / q and a ^ 1, even in 
the framework of the self-similar/stable clustering hypoth- 
esis, where the functions n q (a) are expected to approach 
constants in the highly nonlinear regime (e.g., Davis & Pee- 
bles 1977; Peebles 1980; Balian & Schaeffer 1989, hereafter 
BS). 

In this paper, we use measurements of the PDF by 
Colombi, Bouchet & Hernquist (1996, hereafter CBH) in 
iV-body simulations of flat universes with scale-free initial 
conditions to show that, within the accuracy of the calcula- 
tions, all functions n q superpose to give a unique function 
n ff 



P = PpT(p,eriin 



(4) 



It is important to notice that the measured cumulants in 
the observed galaxy distribution are also in good agreement 
with PT predictions in the weakly nonlinear regime (e.g., 

t Note that we have exploited the scale-free nature of the system 
to express the quantity S q (£, t) as a function of only the variance. 
This would not be possible for scale dependent initial conditions 
such as Cold Dark Matter (CDM). 



nefr(o-) = n q (a). 



(9) 



The measurement of S q is difficult when q is larger than a 
few, because of systematic errors due to finite volume effects 
(e.g., Colombi, Bouchet & Schaeffer 1994, 1995, hereafter 
CBSI and CBSII). However, one can check the validity of 
the property (^) from the overall shape of the PDF rather 
than with its moments. This would mean that PT can be 
applied to the nonlinear regime as follows: 



Ppt[p, er, n e ff(a)]. 



(10) 
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We shall refer to this generalization as Extended Pertur- 
bation Theory (EPT). The only quantities that would be 
needed to compute the PDF, if EPT applies, are the func- 
tions n e fi(a") and the variance itself. If valid, EPT would 
thus result in tremendous simplification. 

This paper is organized as follows. In Section 2, we 
present the measured quantities n q , 3 < q < 5 as func- 
tions of the variance in our iV-body simulations and show 
that they indeed define a unique function n e g(a) that can 
be fitted analytically. We then compare the PDF predicted 
by EPT to that measured for different values of a. The re- 
sults allow an extrapolation to non power-law spectra. In 
section 3 we discuss the asymptotic behavior of the density 
PDFs in the EPT framework which depends only on a few 
parameters. We relate them to the n e s found. We conclude 
in section 4 by summarizing our results and discuss them in 
light of related works. For completeness, in the appendix we 
summarize the results obtained from PT, the generic prop- 
erties expected for the density PDF, and the algorithm used 
to compute it numerically. 



2 MEASUREMENTS 

The scale-free simulations we analyze are described in detail 
in CBH. They were performed with the cosmological tree- 
code of Hernquist, Bouchet & Suto (1991), using 64 3 parti- 
cles and periodic boundaries. We consider four values of the 
spectral index: nii ncar = —2, —1 (two iV-body sets), and 
+ 1. CBH measured the count probability distribution func- 
tion (CPDF) Pn(£) at various times in these simulations. 
The CPDF is just the probability of finding iV particles in 
a spherical cell of radius £ randomly thrown in the data set. 
Thus, the CPDF differs from the PDF: discreteness effects 
have to be accounted for to transform the CPDF into the 
PDF since the latter involves the underlying continuous den- 
sity field. CBH computed the quantities S q , 3 < q < 5 from 
the moments of the CPDF, with the proper corrections for 
discreteness effects. They studied other contamination ef- 
fects, such as artificial correlations related to the way initial 
conditions are generated, short range softening of the forces 
in the iV-body simulations, and finite volume effects aris- 
ing from missing waves longer than the box. In particular, 
CBH defined a reliable scale range over which finite volume 
corrections were applied. From these measurements of the 
functions S q (a), we can directly infer n q {o) for 3 < q < 5. 

This section is organized as follows. First we study the 
functions n q (a) in § 2.1. We will see that they all superpose 
to a single function n c g (a) within the measurement errors. 
In § 2.2, we check whether this result applies to the CPDF 
itself, and in § 2.3, we present a possible extrapolation to 
non power-law spectra. 

2.1 The effective spectral index n c ff 

In Figure |lj we show the dependence of the deduced n q on 
the variance, for each value of nii noar . The various symbols 
correspond to different values of q. They superpose, which 
is not surprising in the weakly nonlinear regime since all 
functions n q should tend to the same value nii near according 
to PT (although not perfectly; this point is discussed in 
detail by CBH). More interestingly, the superposition still 
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Table 1. Parameters used in fit I (eq. Jill). 
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holds in the nonlinear regime where the vertical spread of the 
points does not increase (at least not significantly), which 
corresponds to a very peculiar behavior of the statistics in 
the nonlinear regime, as discussed in the introduction. 

In figure [l], some analytical fits based on two different 
approaches are also shown. In each panel, the thick long 
dashes were calculated from 

nlsia) = 71lm Car + 

(^nonlinear niinear )x T /(x T + x T ), (11) 

with 

x = exp [log 10 (a 2 /cro)] ■ (12) 

Table ^ gives the corresponding values of the parameters 
"-nonlinear, T and o~q as functions of niinear- In equation (|ll[), 
we impose that the index n e s reaches a constant plateau 
at large er, n non ii near . The parameter ero gives the position 
of the transition between the linear value and the nonlinear 
value, and r defines the width of this transition. 

Alternatively it is possible to allow n non iincar to have 
a residual a dependence like the one found in CBH. The 
thick short dashes in Figure |l| use the following interpolation 
formula 

n e s(cr) = n cfi (a) + 

^■nonlinear 

(a)-n l aH (a)]x/(x + l/x). (13) 

The parameter x is given by equation (jl2|) with log 10 o~% = 
1.8, 1.8, 2.2 and 3 respectively for n = —2, — 1, and +1, 
and the function n n0 niincar(o") is 
34 

"nonlinear^) = y - 3 - S 3 (a*/ 100)°° 45 , (14) 

wherej? 3 is the value of S™ cas for a = 10 found by CBH. (We 
have 53 = 10.3, 4.8, 3 and 2.4 respectively for niinear = — 2, 
-1, and +1). 

The two fits (0) and (|l|), designated hereafter respec- 
tively by I and II, are thus based on different theoretical 
prejudices. They both suppose that PT predictions are veri- 
fied in the weakly nonlinear regime a 2 < 1. Model I assumes 
moreover that the stable clustering hypothesis applies for 
a 2 ^> 1, i.e. that the functions S q do not depend on scale 
(or on a) in the highly nonlinear regime. Of course, this im- 
plies that n-eff becomes asymptotically constant at large a, 
and equal to n n0 niinear- Model II takes into account the slight 
deviation from a constant found by CBH for the functions 
gmeas j n ^ ne highly nonlinear regime, namely 

S7= as K <i 2x0 ' 045(? - 2) . (15) 

Note that such a power-law behavior is actually incompati- 
ble with the existence of a universal function n e s{cr). On the 
other hand, this inconsistency is very mild and produces de- 
viations from property (W well below the uncertainties in the 
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Figure 1. The measured spectral indices n q obtained by inverting equations (g) (triangles), (|6|) (squares) and (|ij) (pentagons), using 
measured by CBH. Each panel corresponds to a different value of the initial spectral index nij near . In the case 



the values of 5™ cas 



= —1, two simulations have been used (open and filled symbols). In each panel, the thick long dashed curve uses equation ( |l l[ ) 
with the parameters given in Table The thin short-dashed curves bracket the range of possible fits, given the uncertainties in the 
measurements (see text). The thick short-dashed curves use the prejudice ( jL4| ) for the large a behavior of n non lincar ( see text). 



measurements. Such disagreement could explain, however, 
the slight broadening of the vertical scatter of the points at 
very large a for n = and n — +1. 

Figure [l] also contains two thin short-dashed curves 
above and below the thick long-dashed one. They use the 
same parameterization (|ll|), but with different values of 
nnoniinear, which we denote ra~ oa , Jnear and n+ nlincaT for the 
upper and the lower curves respectively. To some extent, 
they bracket the domain of possible values of n e ff, given the 
uncertainties in the measurements. The numbers n^,. in , jr 

IlOllllllL ci 1 



come from estimates of the errors in the measurements 5'™ cas 
of the parameters S q 
ear regime 



by CBH as follows. In 



the nonlin- 

(a) = 



such errors can be described as S'™ cas ' 
S^. The values of the uncertainty factors f q 
are given by CBH. We take n^^. = (nj) g=3 ,4,5 where 

assuming that S q is cal- 



is the function inferred from S„ 



culated from equations (0) 
with the value of n non ii nea r 



, (^) and (fjj) by replacing ^linear 
given in Table ^. 

The lower panels of Figure |l| contain symbols that 
are significantly outside the range defined by the two thin 
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dashed curves, at a 2 < 1 and for a 2 of the order of few 
units. In the first case, this is because our estimate of the 
errors tends to zero in the weakly nonlinear regime, since 
we assume that n e fi reaches nij near . In the second case, the 
method used by CBH to correct for finite volume effects can 
be used only when a 2 is much larger than unity. As a result, 
it is quite difficult to apply it accurately in the regime where 
a 2 is only a few units. 



2.2 The probability distribution function 

As described in the previous section, we find that the first 
few moments of the PDF are well fitted by the predictions of 
EPT when the spectral index and the variance are allowed 
to be adjustable parameters. To see whether this property 
holds at higher order, we compare in Figure |^ the measured 
CPDF in our A-body sets (tiny squares) to what is predicted 
by EPT (solid curves), for various values of a and average 
number of object per cell 

oo 

(TV) ee J] NP N (16) 

N=0 

(see Table 2). The details of the computations of the PDF 
in the EPT framework are given in detail in Appendix A. To 
obtain the CPDF, one convolves the final result with a Pois- 
son distribution (see, e.g., Coles & Jones 1991). Concerning 
the function n e a(a), we take the fit (|ll| ) with the parameters 
given in Table |l| The dashed curves use the functions n e fr(<r) 
corresponding to the thin short-dashed curves in Figure hi 
To some extent, they reflect the uncertainty in the value of 
n e g. However, there are other significant sources of error, 
which we list here: 

• There is an uncertainty in the value of a used to com- 
pute the CPDF in the EPT framework, since we use a as 
measured in the A-body experiments. For nii ne ar > — 1, 
the error in this estimate is expected to be rather small, 
less than or of order 30%. For nii ncar = — 2, finite volume 
effects strongly contaminate the measurements (CBH). In 
that case, we take the polynomial fit for a 2 (as a function of 
scale) in logarithmic coordinates computed by CBH, which 
is expected to give values of a 2 correct to ~ 30%. We do 
not take explicitly into account the uncertainty in the mea- 
surement of a in what follows, although it is by definition a 
requirement to know a accurately for EPT to give sensible 
predictions. Indeed, since the variance strongly influences 
the shape of the PDF, particularly its width, it is obvious 
that two PDFs with different values of a 2 cannot be ex- 
pected to agree with each other (unless one considers a finite 
range of values of p for the comparison) . 

• There is an uncertainty in the measurement of the 
CPDF itself. In particular, the tails of the CPDF are sub- 
ject to a variety of spurious effects (discussed in detail in 
CBSI, CBSII and CBH): (i) some memory of the pattern 
used to set up the initial conditions is likely to be preserved 
during the simulation, particularly in under-dense regions. 
In our case, particles are initially perturbed from a mesh, 
inducing grid effects which are expected to alter the mea- 
sured CPDF at small values of N, increasing with — niinear 
(CBH). These effects can in principle be avoided by using 



Table 2. Values of a and {N) used for the curves in Figure 2. 



^linear 


curves 


(a) 


(b) 


(c) 


(d) 


(e) 


-2 


(T 


1.49 


5.28 


5.1 


16.2 






(TV) 


1100 


4.37 


0.276 


0.0174 




-1 


(X 


0.57 


0.56 


3.72 


3.36 


17.94 




(N) 


1100 


69.3 


4.37 


0.276 


0.0174 





a 


0.38 


3.67 


3.46 


22.9 






(N) 


1100 


4.37 


0.276 


0.0174 




1 


a 


0.32 


4.19 


4.0 


30.3 






(N) 


1100 


4.37 


0.276 


0.0174 





a "glass" to generate initial conditions^] (White 1994). (ii) 
Because of the finite size of the simulation box, the large- N 
tail of the CPDF is dominated by a few rare clusters, induc- 
ing larger and larger (spurious) fluctuations in the CPDF as 
A increases, followed by a sharp cut-off. 

With this in mind, we see that the agreement between 
EPT and the measurements is excellent, except for nii ncar = 
— 2, a case we discuss further below. There are some very 
small discrepancies at the largest value of a we analyzed 
for nii Ilcar = and 1 [right curves of top panels of Fig. ^, 
with label (d)]. Otherwise, the solid curves almost exactly 
superpose on the small squares. 

For nii ncar = —2, there are some significant disagree- 
ments between EPT and the measurements. While the mea- 
sured high-A tail of the CPDF is quite well described by 
EPT, there is disagreement in the second curve from the 
bottom [label (b)], where the shape of the CPDF predicted 
by EPT appears to be incorrect at small and moderate N. 
For the lowest curve [label (a)], the agreement between EPT 
and measurements is not as bad as it looks, given the un- 
certainties in the function n e g and the large errors in the 
measured CPDF arising from finite volume effects. This is 
actually not surprising since the results are nearly in the 
weakly nonlinear regime where PT is expected to apply. 

To summarize our nii ne ar — — 2 analysis, EPT does not 
seem to describe very well what is happening at small and 
moderate A, i.e. in under-dense regions, but is still quite 
accurate for large values of A, i.e. in over-dense regions. 
This result is not in contradiction with the fact that we 
infer a single function n c ff(er). Indeed, we measure it from 
low-order moments of the CPDF, which are dominated by 
the large- A tail of the CPDF, increasingly so with the order 

q- 

2.3 Application to non power- law spectra 

The previous results apply to scale-free power spectra with 
particular values of niinear. The first step to generalize them 
to arbitrary (Gaussian) initial conditions is to extrapolate 
them to other values of nii noar . Here, we propose a series of 
simple fits for the values of the parameters nonlinear , ao and 

• Of course, this procedure does not correct for the discrete na- 
ture of the simulated data but keeps the white noise level of the 
particles as low as possible, without introducing grid effects (see, 
e.g., Baugh et al. 1995). 
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Figure 2. The measured CPDF as a function of N/(N) (small squares) for various values of (N) and a (see Table 2), where (TV) = 
N NPpf is the average number of objects per cell. Each panel corresponds to a different value of ni; near . The predicted CPDF from 
extended perturbation theory (cq. |Eu|) (solid curve), is calculated assuming that the function n c g (a) is given by equation ( fil| ) with the 
parameters displayed in Table [l]. The two dashed curves bracketing the solid one correspond to the upper and the lower thin dashed 
curves in Figure hi 



t, of fit I as functions of nii ne ar: 

3 (ni incar — 1) 



^nonlinear \J^\ 



3 + m incar 

T(niinoar) = 0.8 - 0.3 n U n 
l0g 1Q (To( n lincar) = 0.2 - 0.1 Ula 



(17) 

(18) 
(19) 



As can been seen from Figure ^[ these fits reproduce well the 
values found for nii ncar = — 2, —1 and 0, but are less accu- 
rate for roii ncar = +1. In any event, the values of cosmologi- 
cal interest are within the range niinear < 0. Figure ^ shows 
the quantities n} ea obtained for different values of the index 
^linear as functions of the variance (solid lines). [The inter- 
polations (|l7|), ( |li| ) and ( |is| ) are used only for ni incar = —0.5 
and nii nC ar = -1.5.] 

Note that Fry, Mellot & Shandarin (1993) measured 
S3 ~ 9/(3 + nii ncar ) in the highly nonlinear regime, which 
implies nonlinear (1.9nu ne ar - 2.5)/(3 + niinear), a formula 
very similar to equation (|17[). This equation diverges for 



wiincar = —3, implying that equation (jl^) is probably valid 
only for -3 < niinear & 0. 

Now, we can generalize the above results to scale depen- 
dent power-spectra. This is actually not completely straight- 
forward. For scale-free initial conditions, a single parameter 
fiiinear can be used to entirely describe the shape of the ini- 
tial power-spectrum. For more complex initial conditions, 
such as CDM, PT predictions require not only a knowledge 
of — (rtiincar + 3) = d log af incar / d log £ , but also of higher log- 
arithmic derivatives of the linear variance with respect to 
the filtering scale 



Iv 



' rr 2 

lineal' 



dlog p £ 



P > 2. 



(20) 



For example, equation (JSJ> must now include a term depend- 
ing on 72 and more generally, the ratio S q involves -y p terms 
forp < 2. It was, however, noticed by Bernardeau (1994b) 
that for CDM like power spectra, such corrections can be ne- 
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Figure 3. The nij near dependence of the parameters n non i; near 
(squares), r (pentagons) and log 10 Op (triangles) used in the fit I 
[eq. (UU)]. The symbols correspond to the values given in Table 2. 
The solid curve, long dashes and short dashes correspond to the 
simple analytical forms (|l^), jlq) and (F9I) respectively. 



glected for S q , q < 4, i.e., one can take y p = for p > 2, while 
still preserving reasonable accuracy. This may not be true 
when higher order moments must be predicted accurately. 
Baugh et al. (1995) took into account 7 P for p < 3 to predict 
S q up to q — 10. But these are tiny refinements and it has 
been shown that, in the weakly nonlinear regime, the shape 
of the PDF is quite well reproduced when only the local in- 
dex nii near is used (B94) . There is no reason why this should 
also be the case in the EPT framework, but we can assume 
that it is so and check this hypothesis directly by measuring 
the functions n q in simulations with scale-dependent initial 
conditions. 

Figure ^ displays the functions n q (a) obtained from two 
CDM simulations with different box sizes (by Davis & Efs- 
tathiou 1988 and Frenk et al. 1990), using measurements of 
the ratios S3 (triangles), S4 (squares) and S& (pentagons) 
by CBSI. Again, as in Figure hi the quantities n q superpose 
to define a unique function n e g(a), in agreement with EPT. 

Now it remains to be seen whether one can infer the 
function n c g (a) from our results for scale- free power-spectra 
[using eqs. (jlT]), @ and @]. The problem is that the 
quantity nii ncar now depends on scale. In the weakly nonlin- 
ear regime, we have of course 

niw(*) = -7i(*)-3, a 2 «l. (21) 

But in the nonlinear regime there is an ambiguity concerning 
the scale at which the initial index should be chosen. We 
consider two simple options: 

(A) With an Eulerian approach, one would generalize for- 
mula to any value of a. 

(B) Within a Lagrangian framework, one would rather 
choose nii ncar according to 



Figure 4. The effective index for the CDM model. The trian- 
gles, squares and pentagons correspond to the measured functions 
n q (a) for q = 3,4 and 5 respectively. The solid lines correspond 
to ourphenomenological fit n^ ff for scale-free initial conditions 
[ eg. (|ll|)] with parameters given by Table |lj and by eqs. ([n]), 
jlq ) and ( fl9| ) for rti; ncar = —0.5 and nii ncar = —1.5. The num- 
bers n on the figure stand for nij noar . The long dashes correspond 
to the prediction for CDM using the prescription (A) and the 
short dashes use the prescription (B) (see text). 

mmeax(^) = -7l(4)-3, (22) 

where the scale £ s is related to I in the same way as in 
the Hamilton et al. (1991) prescription for the transform 
function of the variance, i.e. 

i\ w £ 3 (1 + a 2 ). (23) 

Of course, we have l s — £ when a 2 -C 1, in agreement with 
equation (|2l|). 

In Figure ^, we compare the prescriptions (A) (long dashes) 
and (B) (short dashes) to our measurement of function 
n e ff(o")- The first prescription appears to be clearly better 
than the second. This is somewhat surprising since prescrip- 
tion (B) is expected to be more "physical" than (A). Further 
investigations are required to understand the implications of 
this result. 



3 THE DENSITY PDF IN THE HIGHLY 
NONLINEAR REGIME 

In the continuous limit, and in the highly nonlinear regime, 
one can use PT results to compute some simple analytical 
properties of the PDF in the EPT framework. Indeed, as 
we are going to see in § 3.1, the full shape and properties 
of the PDF in this regime can be simply described by a 
few parameters, which are of course all fixed by the value 
of n e ff- One can thus check EPT predictions by measuring, 
if possible, these parameters directly from the shape of the 



© f996 RAS, MNRAS 000, 



8 Colombi et al. 



PDF, and by comparing the results with what would be 
obtained with the values of n e ff we measured in § 2.1. This 
analysis, which is the subject of § 3.2, should confirm the 
results of 5 2.2. 



3.1 EPT results 

Here, extrapolating the PT calculations of B94 to the highly 
nonlinear regime (i.e. assuming S q = Sq ), we sketch sim- 
ple analytical properties of the density PDF in the EPT 
framework, when the variance is large. Details are given in 
Appendix A. Let us just recall that in the continuous limit, 
the shape of the density PDF is related to two functions 
g(z) and h(x), which are both entirely determined from the 
expression of the generating function of the S q parameters 
(BS). Their respective domains of validity are bounded by 



Po = a 
and 

pu — a 



-1) 



(24) 



(25) 



where the parameters a and u> are set by the value of n B g (see 
Table 3 for numerical values and Appendix A.l for analytical 
expressions of o and u>). There are two regimes: 



P PT (p, a, n cft ) 
and 

Ppt(p, cr, n cff ) 



-9 



pi 



h 



when p < p , 



when p 2> Pu ■ 



(26) 



(27) 



These regimes actually overlap in the region p u <SC p -C Po 
when the variance is large. In this interval, functions g and 
h exhibit the same power law behavior. Note that the global 
shape of g(z) depends only on the value of u> (see BS, B94 
and Appendix A). It exhibits an exponential cut-off at small 
z. The function h(x) is more complex and depends on the 
whole shape of the generating function of S q . It is expected 
to exhibit the following power-law behavior at small x 



h(x) 



r(w) 

and an exponential tail at large x 
3a s 



h(x) 



40F 



5/2 



exp(-\y s \x) 



x> 1. 



(28) 



(29) 



Again, the parameters y s and a s are entirely fixed by the 
value of n e g (see Table ^ and Appendix A.l for the analytical 
expression of y s ). 

Figure |^ shows the functions z 2 q(z) and x 2 h{x) [nu- 
merically obtained from equations (^) and ([I9]) in Ap- 
pendix A.l] in logarithmic coordinates for various values of 
n c ff = —9, —3, —1 and —0.5. These are approximately the 
values of n e g(a) we measured in our scale-free simulations 
at a — 10 for nunear = —2, — 1, and +1 respectively. 



3.2 Comparisons with TV-body Results 

As for the CPDF in § 2.2, we can check the validity of EPT 
by testing directly the analytical forms of § 3.1 in our TV- 
body simulations. In particular, we can try to measure the 



Table 3. Parameters 
spectral index n e ff. 



a, a s and \y s \ for various values of the 



"■eff 




a 


a s 


vi».i 


-0.5 


0.6667 


2.044 


18.94 


1.7436 


-1 


0.6000 


2.104 


6.677 


2.4923 


-2 


0.5000 


2.121 


2.809 


3.9574 


-3 


0.4286 


2.082 


1.842 


5.4112 


-i 


0.3750 


2.025 


1.413 


6.8614 


-5 


0.3333 


1.966 


1.169 


8.3100 


-6 


0.3000 


1.908 


1.012 


9.7576 


-7 


0.2727 


1.855 


0.900 


11.205 


-9 


0.2308 


1.763 


0.752 


14.098 


-11 


0.2000 


1.689 


0.657 


16.991 



parameters fixing the shape of the functions g and h (if they 
exist, of course), and compare the results with what would 
be obtained from the values of n e s we measured in § 2.1. 
This analysis should confirm the results of § 2.2. However, 
it is not easy to carry this out, for several reasons: 

• The effective power index n e g can vary with a even in 
the highly nonlinear regime, which would render the com- 
parisons with EPT quite difficult, since each value of the 
variance should be considered independently. However this 
variation is, as discussed in CBH and in § 2.1, quite weak, 
and n e fl can be taken to be constant to a first order of ap- 
proximation. 

• Strictly speaking, the continuous limit is not reached 
in our TV-body simulations. Discreteness effects can signifi- 
cantly contaminate the measurements, even in samples with 
~ 3 x 10 objects like our TV-body sets. Fortunately one can 
correct for them (at least partially), as proposed by Bouchet, 
Schaeffer & Davis (1991, hereafter BSD). 

• As discussed in § 2.2, other spurious effects, such as 
grid effects at small TV and more importantly finite volume 
effects at large TV can significantly affect the measurements. 

Except for the first point above, these problems have been 
extensively discussed in BSD, Bouchet & Hernquist (1992, 
hereafter BH) and particularly in CBSII. We refer the reader 
to these papers for details concerning the practical measure- 
ment of the functions g and h. 

Although we do not find it necessary to present the 
results here, we find a remarkable qualitative agreement be- 
tween the measurements and the predictions ( prj[ ) and ( p7[ ) in 
the regime a 2 3> 1, after the appropriate corrections. Aside 
from the case nn ne m = +1> and to a lesser extent nii noar = 0, 
the scale dependence of n e g cannot really be neglected (see 
Fig- 0). 

We also tried to estimate by direct measurements (see, 
e.g., BSD, BH, CBSII) the most important parameters of the 
functions g and h, i.e. uj and \y s \. The first one, to, indeed 
completely determines the behavior of the PDF around its 
maximum and in under-dense regions, and the second one, 
\y 3 \, determines the high-density tail of the PDF, therefore 
influencing enormously the values of the moments of the 
PDF. The parameter u) is thus estimated by studying the 
small-TV behavior of the CPDF (through an estimator of 
the function g, using the corrections for discreteness pro- 
posed by BSD), whereas \y s \ is estimated by analyzing the 
large-TV behavior of the CPDF (through an estimator of 
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x 
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Figure 5. The functions z 2 g(z) (left panel) and x 2 h(x) (right panel), in logarithmic coordinates, for various values of n cr j = —0.5, — 1, 
—3 and —9 (solid curves). When n e ff decreases, the exponential tail of the function z 2 g(z) at small z becomes sharper and the power-law 
at large z less steep. The maximum of x 2 h(x) increases with n ef j . The dashed curves are the "preferred" fits of the function h measured 
in our iV-body simulations, using equation (pOh with a, b, c, u> and \y a \ as adjustable parameters. If EPT applies, each dashed curve 
should, within measurement errors, superpose to the globally closest solid curve. 



Table 4. Comparison between EPT and measurements in our iV-body simulations for w and \y s \- 



^linear 


™cff 






y\y s (h)\ 


l/bs(n e fi)| 


-2 
-1 


+ 1 


_o+2.18 
a=c 2.75 
Q+0.67 
O=c 0.79 
i_i_0.32 
i=t 0.36 
n ci 0.24 
u.ozc 2g 


3+0.05 
0.5 ±0.05 
0.65±°;f 

n fic;+0.35 


0.43 ± 0.04 
0.6 ±0.04 
n fi7+0-03 


18±1 
5±o.5 

2± :i 


14±| 
5.4± r 2 
2.5 ±0.5 
1.74 ±0.4 



the function h, using the corresponding BSD corrections for 
discreteness). Table ^| summarizes the results of our mea- 
surements (u(q), l/l?/.,(fe)| ). c ompared with what is obtained 
from equations (|40j) and (|l2|) in Appendix A.l taking for n 
the value of n e g we measure from quantities Sq at a 2 ~ 100 
(w(n e ff),l/|2/ s (n e ff)|). 

The agreement between both estimates of uj and \y s \ 
is quite good, as expected if EPT applies. There is a small 
problem in the case nii ncar = —2, where the measured lo 
seems larger than the one given by EPT (but the error 
bars on each estimate overlap). The same effect exists for 
^linear = — 1, but it is less significant. 

We also tried to find a best fit for the function h, using 
the parameterization proposed by BSD, i.e. 



h flt (x) = o % 7 ^ x^- 2 (l + bx)- c e X p(-\y s \x), (30) 
r(w) 



with a, b and c being adjustable parameters^], and with our 
"preferred" values of ui and \y s \ obtained from the direct 
measurement of the functions g and h. The resulting param- 
eters are listed in Table ^. We obtained the dotted curves 
in the right panel of Figure || which compare well with the 
predictions from EPT, with our "preferred" measured value 
of rieff at a 2 ~ 100, particularly for the large- a; exponential 
tail. The agreement with EPT is not that good at small x 
for n e g = —9 (nii Ilca r = —2) and is marginal for n e g — —3 



§ Note that, in principle, in equation (pc|), the parameters a, b 
and c are entirely determined in the EPT framework. The number 
a is indee d g iven by equation (^j) in Appendix A.l and from 
equation (129), we get 



+ 1/2, 

4 v / 7ra(l - uj) 
3a s r(w) 



1/(^+1/2) 



(31) 
(32) 



However, with such values of a, b and c, the simple analytical 
form (pd) approximates the exact EPT result satisfactorily only 
for n e ff > —2. At the degree of approximation we consider here, 
it is thus better to keep a, b and c as free parameters. 
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Table 5. Our preferred values of a, b, c, uj and \y s \ in equation 
(|30j) as we measure functions g and h. 



^linear 


a 


b 


c 




1/lw.l 


-2 


1.273 


1.455 


0.9 


0.3 


18 


-1 


2.068 


4.3 


0.6 


0.5 


5 





2.7 


8 


0.225 


0.65 


2 


+ 1 


1.937 


0.486 


-1.35 


0.65 


1 



(^linear = — !)• Of course, all these measurements are quite 
delicate and the uncertainties are large. However, the dis- 
agreement for n c ff = — 9 is quite consistent with the conclu- 
sions of § 2.2. 

Thus, the results of this section corroborate those of 
§ 2.2, as expected. 



4 DISCUSSION 

In this paper, we have studied the probability distribution 
function (PDF) of the density field smoothed with top hat 
windows in iV-body simulations with scale-free initial con- 
ditions and linear spectral indices nii ncar = —2, —1, and 
+ 1. We attempted to determine if the following Extended 
Perturbation Theory (EPT) applies. Such an ansatz consists 
of assuming that the PDF is given by P(p) = Ppt(p, cr, n e g), 
where the function Ppt{p, c, nii ncar ) is the prediction of Per- 
turbation Theory (PT, e.g., B94), which is valid in the weakly 
nonlinear regime a ~ (Jii ncar <C 1. The difference with PT 
is that the effective spectral index n c ft is assumed to be an 
adjustable parameter of the model and a 2 is the nonlinear 
variancepj]. The measurement of low-order moments of the 
PDF and of the PDF itself in our iV-body experiments shows 
a remarkable agreement with EPT, except for nii ncar = —2 
in some regimes. (This is however a limiting case, subject 
to many numerical uncertainties, particularly finite volume 
effects, see e.g. CBH). 

As can been seen from Figure jj], the effective spectral in- 
dex, equal of course to nn near in the weakly nonlinear regime, 
decreases with a to approximately reach another plateau at 
large o. If this second plateau were completely flat, our mea- 
surements would agree with the predictions of stable clus- 
tering. This does not appear to be exactly the case, except 
perhaps for nii Ilcar = —2, although the deviation from sta- 
ble clustering is quite weak. We do not discuss this problem 
further here, since it was analyzed extensively in CBH. The 
difference between the second and the first plateau increases 
with —niincar. For a 2 ~ 100, we have for example n c g ~ —9, 
—3, —1 and —0.5 respectively for nii ncar = —2, —1, and 
+ 1. 

It is important to notice that the quantity — (3 + n e fT) is 
from the logarithmic derivative of o 2 with respect 



* Thus, the PDF obtained from EPT has exactly the same vari- 
ance as the data, by definition. This is obviously a necessary con- 
dition for agreement between theory and measurements (but not 
a sufficient one, in general, except of course in the Gaussian limit). 
This is why we think it does not make much sense to try extend- 
ing PT results to the nonlinear regime by identifying the variance 
of the theoretical PDF with the value given by linear theory. 



to scale, except of course in the weakly nonlinear regime. 
Indeed, as measured in iV-body simulations by Efstathiou 
et al. (1988), CBH, and Jain (1996), the variance follows 
quite well the stable clustering predictions in the framework 
of the self-similar solution (see, e.g., Davis & Peebles 1977; 
Peebles 1980), in the regime a 2 J; 100. In particular, we get 



(3 + n CT ) 



dlog a 



3(3 + n lin 



<: 100, (33) 



dlog^ 5 + ni inear 

-1.5, —1.2, —1 respectively for ni; ncar = —2, 
These values are quite different from the n e fr 



i.e., n a = —2, 
-1, and +1. 
we measure. 

We have tried to apply EPT to a scale dependent power- 
spectrum, such as CDM. The results given in Figure ^ are 
encouraging. However, further investigations are needed to 
address the case of non power-law spectra in greater detail. 

Szapudi, Meiksin & Nichol (1996) have recently mea- 
sured the functions n q in the Edinburgh/Durham South- 
ern Galaxy catalogue (EDSGC). In this angular survey, the 
functions n q superpose well, given the uncertainties in the 
measurements, although there is a systematic increase in 
the scatter of values of n q on small scales. Szapudi et al. 
find, as we do, that the function n e fr exhibits two plateaus, 
one on small scales where —7 ^ n e g ^ —4 and another on 
large scales where —2.5 ^ n c ff ^ —1, in good qualitative 
agreement with the results of our analysis and with what 
would be expected in a CDM-like universe. Of course, such 
a comparison is limited by the fact that there are non trivial 
projection effects of the galaxy distribution along the line of 
sight in angular surveys. Moreover, one has to keep in mind 
that light does not necessarily trace mass. 

If EPT is indeed valid, a knowledge of the variance 
a and of the spectral index n a g is sufficient to compute 
the PDF. We have provided analytical fits to the function 
n e s(a). The nonlinear variance o 2 can itself be calculated 
with other nonlinear ansatz, such as the fitting formula 
found by Hamilton et al. (1991) (improved later by Peacock 
& Dodds 1994; Jain, Mo & White 1995; Baugh & Gaztanaga 
1996), provided one generalizes it to the smoothed density 
field. 

Our analysis has been done for flat universes with den- 
sity parameter Q = 1. PT predictions for the ratios S q are 
expected to be only weakly dependent on the value of Q or 
on the cosmological constant A (e.g., Martel & Freudling 
1991, Bouchet et al. 1992, 1995, Bernardeau 1992, 1994b, 
Hivon et al. 1995). But this does not really guaranty the 
existence of a universal function n e g(a) in the case Q 7^ 1 
and/or A / 0. More investigations have to be undertaken in 
this regime. 

Of course, it also remains to explain why the hierarchy 
of cumulants of the PDF looks similar in the weakly nonlin- 
ear and the highly nonlinear regimes. One way to address 
this question would be to study the transition between the 
two regimes, by computing higher order corrections in equa- 
tion Q3p (e.g., Scoccimarro & Frieman 1996; Scoccimarro 
1996). Alternatively, as in the early attempts of Fry (1984) 
and Hamilton (1988), one could study the BBGKY equa- 
tions (see, e.g., Peebles 1980) to try to find clues for the 
determination of the S q parameters in the highly nonlinear 
regime. 
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APPENDIX 

PRACTICAL CALCULATION OF THE PDF IN 
THE FRAMEWORK OF EXTENDED 
PERTURBATION THEORY 

In this appendix, we first study, in the framework of EPT, 
the derivation of the cosmic density PDF, and its properties 
when the variance is large. All the results are discussed in 
term of the index n. For PT, n is simply the initial spectral 
index, n = nii ncar , and, more generally, in EPT, we have 
n = n c ff . In the following we will suppose that n < 0. 

A.l Analytic predictions 

In general the density PDF is given by an inverse Laplace 
transform, 



P(p) dp = dp 



dy 
2-Kio' 2 



exp 



4>{y) , p y 



(34) 



where <f>(y) is the generating function of the S q parameters 

^^(-ir 1 ^ (35) 

(S\ =S 2 = 1). 

In the PT framework and for scale-free initial condi- 
tions, the function is determined by the following set of 
implicit equations^ (see B94): 

y = ~(l + S) ( ' l - 4)/3 [(l + S) 2/3 -l] x 



[(n + 3)(l + g)' 2/a + l-n] 



(36) 



<P = |(l + S) ( "- 4)/3 [(l + e) 2/3 -l] x 

(1 + g) [-n(l + g) 2/3 + (ti-4)] . (37) 

The integral @ is not easy to estimate analytically. In 
general, one must calculate it numerically (see section A. 2). 
The important properties of the PDF, however, depend only 
on the behavior of (fi(y) at large y and near its critical point 
y s (BS). The absolute value of y 3 is the radius of absolute 



II These equations have been obtained assuming that the spheri- 
cal collapse solution is well approximated by Q = (1 — 3<5i/2)~ 3 / 2 , 
where &i is the linear density contrast and Q — 1 is the present 
density contrast. 



convergence of the series ©. From equations © and ©, 
one can deduce that when i/> 1 

<t>{y) ~ ay 1 "", (38) 
with 

-,3/(n-4) 



1 - n 



3(1 -n) 



and 



(39) 



(40) 



On the other hand, the behavior of cj> around y s can be 
written 



4>{y) - 4>s- r s {y - y a ) + a s (y - y a ) 3/2 - 
The calculation of y s yields 

ys = 'lA-^ 2 (A - 1) [(n - l)A - (n + 3)] , 
where 



(41) 
(42) 



A = 



2) 2 - (n-l)(n-4)n(n + 3)] 
(n- l)(n-4) 



1/2 



+- 



(n-l)(n-4)' ^ 

The value of a a is hard to obtain analytically and we com- 
pute it numerically for various values of n (see Table |^) . A 
knowledge of 4> s and r s is not relevant here. 

We can now define, following the formalism of BS, the 
typical density of an under-dense region as 

p u = a " >a /v 1 , (44) 

and the typical density of an over-dense region as 

Po = o 2 . (45) 

For regions with small enough density, the PDF can be 
written as 



Pp T (p, a, n) ~ — g ( — 
Pu V pu 



P < Po- 



(46) 



The shape of the function g(z) depends only on the value of 
the parameter ui (BS) through the following integral 



9{z) 



duexp[— zu + it 1 " cos 7rw] x 



(47) 



It has a exponential tail at small z and a power-law behavior 
in z"~ 2 at large z (see B94 or BS for details). 

For regions with large enough density, we have 

Pp T (p,(7,n)~ -jftf ) , p^>pu, (48) 

po \po J 

where the function h(x) is related to function cf>(y) through 
the following transform (BS) 

1 f + °° 

h( x ) = ~— du <j>{iu) e mx '. (49) 

27T / 

It exhibits a power-law behavior (^8|) at small x (this regimes 
overlaps the large z regime of function g) and an exponen- 
tial tail at large x given by equation (E9). Note that the 
parameters S q are simply the moments of the function h, 
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Figure 6. Example of paths in the complex plane for y. The 
parameters correspond to n = — 3 and a = 1. The paths cannot 
cross the real axis for y < y s . When the density is small, the 
paths are on the right side. When the density is large, they are 
on the left side and cross the real axis at y = y s . 



x q h(x)dx, (50) 



and we have the following normalizations 



xh(x)dx = J x h{x)dx = 1. (51) 
Jo 



A. 2 Numerical Integration 

The starting point of the numerical calculation of the density 
PDF is equation ( |34[ ) . Since the intervening quantities are all 
known as functions of Q, this latter is our natural variable of 
integration. The expression for the density PDF thus reads, 



P(p)dp 



dp 



exp 



dQ dy(g) 



2nia 2 

4>{y) p_y 

a 2 a 2 



(52) 



The difficult part of the integration is to choose the proper 
path in the complex plane. The original path for y runs 
from — ioo to +ioo along the imaginary axis. But as the 
functions r(y) or 4>(y) are not analytic in y, there is at least 
one singularity on the real axis for y = y s (< 0), and it is 
impossible to move the path to the left of this singularity. 

Actually, the starting point of the path for the numer- 
ical integration is the position of the saddle point, j/ sa ddio 
defined by 



= P - 



d0(y) 



dy 



(53) 



path is chosen in such a way that <f>(y) + py is kept real 
to avoid unnecessary oscillation of the integrand. Using this 
prescription, followed with an adaptive step, the paths that 
are found are similar to the ones shown in Figure m. 
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